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Abstract 

The relative phase between two uncoupled BE condensates tends to attain 
a specific value when the phase is measured. This can be done by observing 
their decay products in interference. We discuss exactly solvable models for 
this process in cases where competing observation channels drive the phases 
to different sets of values. We treat the case of two modes which both emit 
into the input ports of two beam splitters, and of a linear or circular chain 
of modes. In these latter cases, the transitivity of relative phase becomes an 
issue. 
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I. INTRODUCTION 



Since the first observation of Bose-Einstein condensation, the formation and the nature 
of the relative phase between two condensates has been a central issue of many theoretical 
and experimental studies. It has been predicted by Javanainen and Yoo [1] and observed 
by Andrews et al [2] that two interfering Bose condensates exhibit a clear spatial interfer- 
ence pattern. This shows that in a single run of an interference experiment, they manifest 
themselves as being coherent. Furthermore, it was predicted in [1] that two cases should be 
distinguished. When a cold cloud of atoms is first split into two modes, which are separately 
cooled further into two condensates ( "cut - then- cool"), two independent condensates arise. 
Alternatively, two correlated condensates arise when a single condensate is split into two 
parts (" cool - then - cut") [2,3]. Interference pattern from two independent condensates 
can be different for each realization of interference experiment, while correlated condensates 
show the same interference pattern for each run. Cirac et al [4] showed by analytical ar- 
guments that a system consisting of two independent Bose condensates evolves into a state 
with a fixed relative phase if one detects the emitted bosonic atoms while observing their 
spatial interference pattern. 

A number of authors have studied the possible manipulation of phase coherence and 
entanglement between two or more Bose condensates, with tunneling interaction as the key 
mechanism [5-7]. A scheme has been proposed to use an interferometric scheme including an 
atomic beam splitter to recombine two modes in order to reconstruct the state of a two-mode 
condensate [8]. The buildup of a relative phase between two independent condensates has 
also been investigated in the situation that the atoms emitted from the two condensates are 
mixed in a 50:50 beam splitter [9,10]. Two initially independent bosonic modes, described 
by a factorized state, evolve into an entangled state of the two modes after a large number 
of detections in the output ports of the beam splitter. The relative phase distribution shows 
two narrow peaks, at positions determined by the settings of the beam splitter. The most 
probable detection history has the form of bosons bunching into a single output channel. 
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An exactly solvable analytical model has been discussed [10], which allows one to get closed 
expression for the particle detection statistics over two output channels of the beam splitter 
for a fixed total number of detections. It is remarkable that even though both detection 
channels are identical, in the most probable history all particles are detected in the same 
port. This is obviously connected to the bosonic nature of the particles, for which boson 
accumulation applies. This can likewise be interpreted as a spontaneous selection of a single 
relative phase. When the first particle chooses randomly one of the two output ports, the 
following particles have a tendency to choose the same port, and the relative phase of the 
modes converges to one of the phases imposed by the beam splitter. This can also be viewed 
as an example of spontaneous symmetry breaking [11]. The role of interparticle interaction 
is also discussed, and it has been shown that it leads to collapse and revival of the relative 
phase distribution, thereby reflecting the discrete nature of the states of the system [9]. 

In the presence of a single beam splitter, the relative phase converges eventually to a 
single value. It is interesting to consider cases where more detection channels are present 
which tend to project the relative phase on different values, so that a detection from one beam 
splitter favors phase values that are incompatible with the setting of another one. In the 
present paper we consider a number of model cases where such a conflicting tendency arises. 
This raises the question whether in the end the system simply settles down in one of the 
possible phase values, or whether it continues to shift between values, without ever coming to 
a final decision. We consider cases where the detection statistics can be solved analytically. 
Also we study the effect of a direct Hamiltonian coupling between the condensates on both 
the detection statistics and the corresponding behavior of the relative phase. Examples 
of such couplings are tunneling between condensates in two spatially separated potential 
wells, or stimulated Raman coupling between two condensates corresponding to two different 
internal states [12]. We treat the condensates just as modes of bosonic particles, so that 
most of the considerations hold just as well for photons in cavities. 
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II. QUANTUM STATES OF TWO BOSON MODES 



It will be convenient to express the states of two boson modes in terms of spin-coherent 
states (SCS), which is normally defined for the 2 J + 1-dimensional manifold of states with 
angular momentum J [13]. The spin-coherent state \9,<p) is the eigenstate of the com- 
ponent ~vt ■ ~S of the angular momentum vector with the maximal eigenvalue J, where 
it = x cos (j) sin 9 + y sin <p sin 9 + z cos 9 is the unit vector in the direction specified by the 
spherical angles 9 and 0. This state is obtained from the eigenstate of J z with eigenvalue 
J after performing the appropriate rotation. In the context of two boson modes (or two 
harmonic oscillators), an SU{2) representation arises by introducing the fictitious angular- 
momentum operators 

J x = \ (a f 6+?a) ,Jy = Yi ftb-Va),J z = l - (a^a-Vt) , (1) 

where a and b are the annihilation operators for modes A and B. This is the well-known 

Schwinger representation. These operators obey the standard commutation rules of angular 

momentum ([J x , J y ] — iJ z , etc.), so that the matrix form of the operators (1) on the eigen- 

vectors of J z and J attains the shape that is well-known from angular-momentum algebra. 

. — a ^ _ ^ ^ 

Notice that 7 = y(y + 1), with iV = a^a + tfb the number operator. The eigenvectors of J z 

and J are just the double Fock states \n a ,nb). A given number of particles N corresponds 

to the value J = N/2. The eigenstate of J z with this same eigenvalue is the Fock state 

| AT, 0) , so that the SCS with direction ifcan be defined by the rotation 

\e,<f>) N = R(e,<f>)\N,o), (2) 

with the rotation operator 

R(9,4>) = exp(—i(f)J z )exp(—i9Jy)exp(i(j)J z ) = exp[— i9(J y cos0 — sin </>)]. (3) 

The SCS can be represented as a point on a sphere of radius J, specified by the polar angle 
9, and the azimuthal angle <fi. This sphere generalizes the Bloch sphere, describing the state 
of a spin 1/2, or the Poincare sphere which describes the polarization state of a light beam 
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or a photon. In the present case, the radius specifies the number of particles N = 2J. An 
explicit expansion of the SCS (2) in the Fock states follows then from the transformation of 
the creation operators 

R(6, 4>)a ] R\6, 0) = a f cos °- + # sin = ${B, 0). (4) 

The SCS (2) is found after operating N times with the operator c^(9,<p) on the vacuum 
state, which leads to the explicit result 



N I »r \ 



n=0 



N 



cos n - sin""" - e i{N - n)4, Ira, N-n). (5) 
2 2 v 7 



This demonstrates that the SCS |6>, <j)) N can be viewed as a number state in the mode that 
is a linear combination of the modes A and B, and for which the operator 0), defined 
in (4), is the creation operator. In the SCS, the distribution of the iV particles over the 
two modes is binomial, and the angle 6 specifies the average partition by (n a ) = N cos 2 f 
and {rib) — N sin 2 |. The azimuthal angle represents the relative phase. This quantity is 
complementary to the number difference a^a — b%. Number states with all particles in the 
mode A are represented by the North pole of the Bloch sphere (6 = 0), while the South pole 
represents the SCS with all iV particles in mode B. Points on the equator (6 = 7r/2) stand 
for states with equal population of the modes. Since the state (2) (or (5)) is eigenstate of 
N, the absolute phase is fully undetermined. 

The relation between the SCS and the more common Glauber coherent states (GCS) is 
easily found by representing the latter ones in the form 

\r a e-^ a , r b e-^ b ) = e - (r ' +r ^ )/2 ^ J^r.e"^ + ne'^f \vac) . (6) 

N 

These states are eigenstates of a and 6, and they are obviously factorized, so that they carry 
no entanglement between the modes. It is easy to check that they are related to the SCS 
by the expansion 

\r a e-^,r b e-^) = e"^ £ ^ ^ e -^ ^ ^ > (7) 

N v N . 



with the parameters R, 9 and determined by R 2 = r^ + r 2 , tan | = rb/r a , and = a — 0t- 
This indicates that the GCS has a Poissonian distribution of the total particle number N, 
with average value (N) = R 2 , while the absolute phases a and 4>b of both modes are well- 
specified. For bosonic atoms, states with different total number of particles do not superpose, 
according to the superselection rule, so that we have to restrict ourselves to density matrices 
that are diagonal in N . Since the particle number is conjugate to the overall phase, we 
introduce the density matrix 

p(R, M) = ^ fj d<t>a |r a e-^, r 6 e- l( ^-^) <r a e"^, ne'^'® | (8) 

as the uniform mixture of the GCS (6) over the overall phase a , for a given value of the 
relative phase = a — Applying eq. (7) leads to an expansion of this same density 
matrix in the SCS, in the form 

p(r, e, 0) = e~ R2 £ j^r 2N \e, 4>) N N (e, <p\ . (9) 

N 

The density matrix p(R, 0, 0) is therefore diagonal in the particle number N. 

We observe that to each pair of spherical angles 9 and 0, or, equivalently, to each real 
Cartesian unit vector 1?, corresponds a density matrix p(R, 9, 0), and an annihilation opera- 
tor c{9, 0), as defined in (4). Now consider the annihilation operator c(6*o, 0o), corresponding 
to the unit vector T? - 

In this paper we shall use density matrices that can be represented as a superposition of 
the states (9) for a single value of the strength parameter R, in the form 

J <mf{9A)p{Ri8A), (io) 

where we use the abbreviation J dVt = J Q 27r d<p d9 sin 9 for the integration over the Bloch 
sphere. When we express p(R,9,(f)) as in eq. (8), it becomes clear that eq. (18) is just the 
two-mode version of the Glauber-Sudarshan diagonal coherent-state representation of the 
initial density matrix [15], where the P-distribution is uniform in 0a, and is non-zero only 
for a single value of R. This state is normalized as soon as the distribution / is, which we 
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shall assume. Another special case arises when the function / is nonzero only for a single 
value of 9, and uniform in 0. Then the density matrix (10) can be written as 



It follows from the coherent-state representation (8) that in this case the density matrix 
factorizes into a product of separate density matrices for the two modes, implying that the 
state (11) not entangled. The phase of both modes is uniformly distributed, and the state 
is diagonal in both particle numbers n a and n^. 

III. DECAY AND DETECTION STATISTICS OF TWO BOSON MODES 

A. Master equation and detection histories 

We assume that particles are leaking out of the two boson modes A and B, at a total 
loss rate T. The emitted particles are detected after passing through a beam splitter. For 
simplicity, we assume perfect detection efficiency Moreover, their evolution is governed 
by a Hamiltonian H that is supposed to commute with the number operator N, and which 
describes the energy per particle, and possibly interparticle interaction or tunneling between 
the modes. Since the two modes form an open system, their evolution can be described by 
a quantum master equation [14,15] for the two-mode density matrix p, which we formally 
express as 



Here Cq describes the coherent evolution of the system, which is determined by the Hamil- 
tonian evolution, and the loss of the probability of states due to the emission of particles. 
Its explicit form is given by its action on a density matrix 




(11) 




(12) 
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(13) 



while the compensating probability gain is accounted for by 
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Cxp = Y (apa ] + bpY (14) 



For simplicity the loss rate of the two modes is taken to be the same. The solution of (12) 
describes the evolution of the system averaged over all possible detection histories. In fact, 
we are interested in the conditional evolution for specific histories, where the arrival times for 
particles at each detector are specified. Depending on the specific setup, we have to separate 
the total gain term (14) in terms corresponding to each detector separately, in accordance 
with the method of quantum trajectories [9,4,10]. For instance, when a detector is directly 
coupled to each mode, the term a'pa) describes the effect of a detection of a particle from 
mode A, which corresponds to the annihilation of a particle from this mode. Now we consider 
the setup sketched in Fig. 1, where each mode emits particles into the input port of two 
different beam splitters. Detections in the two output ports of beam splitter I correspond to 
the detection operators c± = (a±b)/y/2, and detections in the output ports of beam splitter 
II correspond to the detection operators d± = (a±e~^b)/y/2. The relative phases can be set 
either by using dephasers, or by differences in the pathlengths of the channels. Notice that 
the detection operators are annihilation operators corresponding to a spin-coherent state 
that is represented by points on the equator of the Bloch sphere. For this setup the gain 
operator d can be separated into four terms corresponding to the four detectors as 
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dp = — (c+pc j + + c_pcl + d + p$ + + d-pdVj = 2^2 ^° c s = S ^ ls P- ( 15 ) 

The integral form of the master equation (12) 

p(T) = e CoT p (0) + V f dt e Co ^diP(t) (16) 

allows us after iteration to express the density matrix as a summation and integration 
over detection histories. The contribution to p(T) from the history with detections at the 
successive time instants t\ < t 2 < . . . < th by the detectors s±, S2, • • ■ sl in the time interval 
[0, T] is described by the operator 

p L ({U, s t } , T) = e ^T-t L . 1 ) CuLe c (t L -t L - 1 ) . . . C lsi e c ^p(0) . (17) 



The effect of the detection operators Cu is a sudden change in the density matrix, which 
indicates the quantum-jump nature of a detection. 



B. Detection statistics and phase distribution 

As initial state p(0) of the system we take a density matrix of the form (10), so that 

p(0) = J dnf(9,<f>)p(R,9,<f>). (18) 

When the Hamiltonian only attributes a fixed energy per particle, its form is H = hu)N. 
Since all density matrices that we shall encounter are diagonal in the total number of par- 
ticles, the Hamiltonian has no effect, and can be ignored. The coherent evolution of the 
density matrix is easily obtained from the identity C \<f>, 9) NN (9, <f>\ = —TN \<f>, 9) NN {9, 0|, 
which when substituted into eq. (9) gives the result 

e CoT p(R, 9, 0) = exp[-i? 2 (l - e - rT )}p(R e - rT/2 , 9, 0). (19) 

This shows that the evolution of the density matrix during a detection-free period of time 
only gives a damping of the strength parameter R, without changing the distribution over 
the Bloch sphere. The action of the detection operators on the density matrix is most easily 
obtained by using eq. (8). The action of the annihilation operators on the SCS is found to 
be given by 

a \9, <f>) N = VNcos 6 - \9, ,b\9, <P) N = v^sin °-e^ \9, < P) N _ 1 (20) 

Then a direct calculation shows that 

c(9 , <p )p(R, 9, <p)c\9 , O ) = l -R 2 {\ + vt ■ Tt )p(R, 9, 0), (21) 

with c defined in eq. (4). The unit vectors 1? and Tt in eq. (21) are defined to point 
in the directions specified by the angles (9, 4>) and (9 , <f> ) respectively. This indicates that 
for these operators cpff is proportional to p. The proportionality factor takes the maximal 
value R 2 when the two directions ~Tto and it coincide, and it is zero when the directions 
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are opposite. It is not surprising that this factor depends only on the inner product of the 
two unit vectors, and thereby on the distance between the two points on the unit sphere. 
Application of (21) leads to the expression 

£ ls p(R, 9, 0) = TR 2 g s (9, <f>)p(R, 9, 0), (22) 

where the functions gi for the detectors 1 and 2 are given by 

1 1 

gi (9,(p) = -(1 + sin 9 cos 0), g 2 {0,<f>) = -(1 -sin#cos0), (23) 

and for the detectors 3 and 4 by 

93(0, 4>) = ^(l + sin0cos(0-O),<?4(0,0) = ^(1 -sin0cos(0-O)- (24) 

The functions are determined by the inner product of the unit vector ~Tt, indicated by 9 and 
0, and the unit vectors Tt corresponding to the detection operators c s . These four unit 
vectors are all defined by #0 = n/2, whereas O = and n for s = 1 and 2, and O = £ and 
£ + 7i for s = 3 and 4. The functions g s add up to 1, so that the total gain operator L\ 
when acting on p(R, 9, 0) just gives the factor TR 2 , as it should. According to eq. (22), the 
effect of the ith detection at time U by detector is that the distribution over the Bloch 
sphere is multiplied by the factor g s ., while an overall factor TR 2 exp(— Tti) has to be added. 
In brief, the detection-free periods produce a damping of R, and the detection modify the 
distribution over the Bloch sphere by a multiplication with a function g s .. For a given value 
of the ratio (n a ) / (rib), as specified by the angle 9, the factors g s modify the distribution 
over the relative phase 0, with a contrast that is maximal when both modes contain the 
same number of particles (6 = ir/2). 

The eqs. (19)-(24) allow one to evaluate explicitly the density matrix (17) corresponding 
to a given detection history, with the initial state determined by (18). The contribution (17) 
to the density matrix is then found as 

L 

p L ({U, Sl } , T) = exp[-i? 2 (l - e" rT )] H(TR 2 e-^) 
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x / <mf{9A) 



p{Re- VT '\9A), (25) 



_s=l 

with n s the total number of detections in channel s (with J2n s = L). This contribution (25) 
does not depend on the specific order of the detections in the various channels. The trace 
of (25) specifies the probability distribution of the detection history {tj, Sj} in the factorized 
form 

L 

w L ({t h s t } , T) = F({n s }) exp[-i? 2 (l - e - TT )] ]J(rR 2 e- Tti ), (26) 

i=i 

with 

F({n.})= I cMf(9, 0) fl g? (9, 0) (27) 

J s=l 

the probability that L successive detections occur in the specific order (si, S2,- • -,sl)- This 
factor F only depends on the number of detections n s for each channel, not on the time 
ordering of the detections. The remaining time-dependent factor in (26) is the probability 
density for detections at the specified instants of time, irrespective of the detection chan- 
nel. The conditional density of the system, given the detection history {tj,Sj}, is equal to 
Pl ({U, Si} ,T) /wl ({U, Si} ,T), which is the normalized version of (25). From the expres- 
sion (26) of the probability density one obtains the probability p({n s },T) that in the time 
interval [0, T] there were n s detections in channel s, (s = 1,. . .,4), irrespective of the order 
of the detections. This requires an integration over the ordered detection times, and a mul- 
tiplication with the number of possible orderings of the L detections over the four detectors, 
given the partition {n s }. The result can be expressed as 

p({n a },T)=P L (T)p L ({n a }), (28) 

where Pl(T) gives the probability that precisely L detections occurred in the time interval 
[0, T], irrespective of the detection channel. This distribution is Poissonian with average 
R 2 (l — e~ TT ). The factor PL{{ n s}) is the probability that the L detections are distributed 
over the four detectors by the partition {n s }, and takes the form 
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L! 



(29) 



This distribution is independent of the strength factor R, the detection time T and the decay- 
rate T. Notice that both the distribution Pl{T) over the total number L of detections, and 
the distribution p L ({n s }) of the L detections over the partitions are normalized. 

In summary, we notice that the decay process only has the effect that the strength factor 
R is damped. The effect of a detection is that the distribution over the Bloch sphere is mul- 
tiplied by one of the factors g s , which changes both the distribution over the relative phase 
and the probability distribution for subsequent detections. The probability distribution of 
L detections over the four detection channels is given by (29). After a detection series given 
by the partition {n s }, the normalized distribution function over the Bloch sphere is given 



When the detections in the channels 3 and 4 are ignored, and M detections have occurred 
in the channels 1 and 2, the distribution of these detections over the two channels can be 
evaluated in the same fashion. The result is 



with rii + n 2 = M. The factor 2 M is needed to ensure normalization, since gi + g 2 — 1/2 in 
this case. This expression is a simple generalization of the result of [10] for the case of two 
decaying modes observed through a single beam splitter. The generalization consists in the 
fact that the populations of the two modes need not be the same in eq. (30). Intuitively 
it is obvious that the partial statistics of detections in the channels 1 and 2 is not affected 
when for some reason the detections in the channels 3 and 4 are simply added without 
distinguishing them. This situation is equivalent to the case that beam splitter II is missing, 
and a single detector is just collecting particles in both of its input channels. For a total 
number L of detections, the probability of having rt\ and n 2 detections in channels 1 and 2, 



by the f(9,<P)Yl s gs°(d,<P)/ F({n s }). 



C. Sum rules 




(30) 
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with n\ + n 2 = M < L can be expressed as a marginal distribution of pL({n s }) in which the 
sum n 3 + n 4 is fixed. When using that g 3 + g± = 1/2 we find 

E Pd{n s }) = ^2-VK,n 2 ), (31) 

n 3 +ri4=L—M ^ ' 

with pm the distribution (30) over channels 1 and 2, regardless the detections in channels 
3 and 4. Ths confirms that the relative distribution of the detections over the first two 
channels remains unaffected by the detections in the channels 3 and 4, provided that these 
are not distinguished. 



D. Special cases 

We have noticed that the effect of detections on the phase distribution is strongest when 
the average number of particles is the same in both modes, so we consider the case that the 
polar angle is 6 = 7r/2, so that r a = r& = R/V~2 = r. For this situation, the two-channel 
distribution (30) has been evaluated in ref. [10]. When the relative phase <f> has a well-defined 
value O , the two-channel distribution is binomial 

PM(n 1 ,n 2 ) = (^cos 2 -^sin 2 -^, (32) 

where the most probable detection history has the values n\ = Mcos 2 (0o/2), n 2 = 
Msin 2 (0 o /2). When the phase distribution is uniform, the two-channel distribution was 
found as [10] 

pM(ni ' n2)= ^(lT)( 2 r 2 2 )' (33) 

which displays boson accumulation, with the most probable history specified by (n 1? n 2 ) = 
(M, 0) or (0,M). After such a history, the relative-phase distribution is proportional to 
cos 2M | or sin 2M |, which peaks at the positions corresponding to the output channels of 
the beam splitter /. 

Now we turn to the detection statistics over the four channels when the initial density 
matrix is specified by eq. (11), with equal population of the two modes, and initial uniform 
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relative phase. Then the initial density matrix (18) is equivalent to the factorized form 
p(0) =p a ® Pb, with 

P- = ^J d( t>» \ re ^ a ) ( re ^ a \ > ( 34 ) 

and a similar expression for p b . Both modes have a density matrix that is diagonal in 
the number state, with a Poissonian distribution. Intuitively one would expect that both 
two-channel distributions (32) and (33) are contained in the margins of the four-channel 
distribution PL({n s }), which must be equal to the product of the marginal distribution 
(31) with M = rii + n 2 , and the conditional distribution Pl(ti3, n^ni, 712). We look for 
detection histories with maximum probability. First we notice that the emission proba- 
bility onto both beam splitters I and II is the same, so that for a total of L detections 
a most probable history must have n\ + n 2 = n 3 + n 4 = L/2. (We assume that L is 
even for simplicity.) If nothing is specified on the distribution of the L/2 detections in 
the channels 3 and 4, the distribution over the two channels 1 and 2 is given by eq. (33) 
with M = L/2, with the most probable partitions (ni,n 2 ) = (L/2, 0) or (0, L/2). The 
relative phase has then converged to the value = or = n, which makes the distribu- 
tion over the L/2 detections in channels 3 and 4 binomial. For example, for the partition 
(711,712) = (L/2,0), the partition over the two other detectors has maximal probability for 
(713,714) = (L/2)(cos 2 (£/2), sin 2 (£/2)). Since the pair of detectors 1 and 2 is fully equiv- 
alent to the pair 3 and 4, another history with the same maximal probability occurs for 
the partition (n 3 ,n 4 ) = (L/2,0), with (711,712) = (L/2)(cos 2 (£/2), sin 2 (£/2)). This corre- 
sponds to a relative phase converging to the value — £. In summary, we expect four 
most probable histories for L detections. The partitions over the four detectors attain 
the values (m, n 2 , n 3 , n 4 ) = (L/2)(l, 0, cos 2 (£/2), sin 2 (£/2)), (L/2)(0, 1, sin 2 (£/2), cos 2 (£/2)), 
(L/2)(cos 2 (£/2),sin 2 (£/2),l,0) and (L/2)(sin 2 (£/2), cos 2 (£/2), 0, 1), while the phase has 
converged in these cases to the values = 0, ir, £ and £ + ix, respectively. These consid- 
erations are backed up by a numerical calculation of the probability distribution Pl({^s})> 
for L = 40, equal population of the two wells (6 = 7r/2), uniform distribution over the rel- 
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ative phase 0, while the setting of the two beam splitters is maximally different (£ = n/2). 
The distribution for equal number of detections through both beam splitters is plotted in 
Fig. 2. The most probable histories are marked. The gradual transition between the two 
distributions (32) and (33) is noticed along the axis ni, when varies from (binomial 
distribution over n\ and n 2 = L/2 — ni) and L/2 (accumulation distribution (33)). 

IV. DETECTION STATISTICS OF TWO COUPLED BOSON MODES 

A. Pulsed coupling between modes 

In this secton, we consider the case that the particles emitted by the two boson modes 
A and B are detected directly, without the use of beam splitters, as sketched in Fig. 3. 
Therefore we separate the gain operator in the master equation (12) as L\ = £± a + Cu, 
corresponding to the two terms in (14). The coherent-evolution operator £ is given by eq. 
(13), where the Hamiltonian H describes coupling between the two modes by tunneling, in 
the form 

H = -— (a^b + ab^ = -h5J x . (35) 

In realistic cases we can imagine that the coupling can be switched on during a time interval 
r, which is sufficiently small so that decay during the coupling is negligible. This means 
that the initial state for the decay process is found by applying the pulse evolution operator 

Uq = exp(—iHr/h) = exp(i5rJ x ). (36) 

In the picture of the Bloch sphere, this is a rotation about the x-axis in a negative direction 
over an angle 5r. When the initial state before the coupling is given by (10), the state after 
switching-off the coupling at the beginning of the detection period is 

p(0) = J cMf(9, <f>)U p(R, 9, <j>)Ul (37) 

The contribution to the density matrix from a given detection history {tj, s^} is expressed by 
eq. (17), where now the indices s of the jump operators C\ s can take the values a or b, and 
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where eq. (37) specifies the initial density matrix. The evolution during the detection-free 
periods is given in eq. (19). The effect of the jump operators on the rotated density matrix 
can be expressed using the identity 

C la U Q pUl = YU c a pc\Ul, 

and a similar expression for C\b, where we introduced the counterrotated operators c a = 
UqClU and c b = £/qo£/o- Their explicit expressions are then 

^ ^ St -~- 5r ^ ^ St ^ 5t 

c a = a cos h 20 sin — , c b = msm h ocos — . 

2 2 2 2 

They correspond in the sense of eq. (4) to the two unit vectors Tt a = — y sin St + 2; cos 5t, 
Jt b = y sin 5t — z~ cos 5t, which arise when the opposite rotation is applied to ±z. By using 
eq. (21), the action of the jump operators C\ a and C\ b in a detection history is given by the 
relation 

C la U o p(R,0,<f>)Ui = rR 2 g a (e,<f>)U p(R,e,<f>)Ul 
C lh U p(R, 6, <j>)ut = TR 2 g b (6, </>)U p{R, 9, 01 (38) 

with 

9a{0A) = ^(l+*?-*?a),9b(0,<i>) = ^(l+tt-ttb). (39) 

Notice that these factors add up to TR 2 . The contribution to the density matrix arising 
from the history {tj, Sj} is now easily found in the form 

L 

p L ({U, s,} , T) = exp[-it> 2 (l - e~ TT )] J](ri2 2 e - rt ') 

i=i 



x 



J cMf(9, <f>)g n a "(9, 4>)g n b b {e, <f>)U p(Re- rT/ \ 0, fiul (40) 



which looks quite similar as eq. (25). The probability distribution for detection histories is 
given by the trace of (40), and the detection statistics can be obtained in the same way as 
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above. In analogy to eq. (28), the probability p(n a ,n b ,T) that in the time interval [0, T] 
there were n a detections in channel a, and n b in channel b, irrespective of their order, is now 



where, as before, Pl(T) is the Poissonian distribution of the total number L = n a + n b of 
detections in the interval [0, T]. The factor pL(n a ,n b ), which represents the probability that 
the L detections are partitioned over the two detectors as (n a ,n b ), is 



As an example, we consider the case that before the coupling period the two modes are 
fully decoupled, with equal population, so that the function / is uniform over the equator of 
the sphere. The density matrix before coupling has then the form (11), with 6 = ir/2. When 
moreover the pulse duration is chosen such that 5r = n/2, we find lt a = —y, Tt b = y, and 
the functions g a and #6 at the equator are found as g a ((f>) = (1— sin0)/2, gb{4>) = (l+sin0)/2. 
The distribution pi(n a , n b ) is now exactly the same as in the case of an initally uniform phase 
distribution, with detectors are placed in the output channel of a single 50% — 50% beam 
splitter [10], and we recover the bunching distribution 



The most probable history of L detections is (n a , n b ) = (L, 0) or (0, L). This is understand- 
able, since the relative geometry on the Bloch sphere of the initial state and the detection 
operators is the same in both cases. In the case of detections through the beam splitter, 
the distribution / is initially uniform over the equator, and the detectors correspond in the 
sense of eq. (4) to two opposite points on the equator. A typical detection history then 
projects the phase distribution onto a narrow peak located at either one of the detector unit 



p(n a ,n b T) = P L (T)p L (n a ,n b ), 




(41) 



with 




(42) 
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vectors. In the case of detectors attached directly to the two modes, the mode coupling 
prior to the detections rotates the uniform distribution over the equator about the x-axis 
over an angle n/2, so that the initial distribution before the detection series is uniform over 
the large circle in the xz-plane. The detectors represented by the detection operators a and 
b correspond to the unit vectors ±z, which again are two point on opposite sides of the large 
circle representing the initial state. However, the physical situation is quite different in the 
two cases. For the initially uncoupled states and detection through the beam splitter, the 
relative phase distribution starts out uniform, and it is converted into a narrow distribu- 
tion during a typical detection history in the output channels of the beam splitter. For the 
initially coupled modes and the detections without the beam splitter, the relative phase is 
initially rather well-determined around = and = ir. A typical detection series now 
projects the state of the system onto the state with all particles either in mode A or in 
mode B, with an undetermined relative phase. If at the end of the detection series a second 
pulsed coupling is applied as described by the operator U , the final state after this pulse 
has a well-determined relative phase. The net result of the entire scheme of pulsed coupling, 
detection series and second pulse is the same as the result of just a detection series through 
the beam splitter. In this sense, the pulsed coupling can be viewed as a replacement of the 
beam splitter. 



B. Continuous coupling between modes 

The situation is different when the coupling between the modes is present continuously. 
Then in expression (13) for the coherent-evolution operator, the Hamiltonian is given by eq. 
(35). Since the Hamiltonian commutes with the number operator N, the decay terms are 
not affected the Hamiltonian evolution, and eq. (19) is replaced by the modified form 

e CoT p(R, 9, 0) = expf-i^l - e- VT )]U{T)p{Re- VT ' 2 , 9, 4>)fP(T) (43) 

with U(T) = exp(—iHT/h) = exp(i5TJ x ). The effect of the Hamiltonian on the density 
matrix for a detection history {tj, Sj} can be expressed in the Heisenberg picture, with the 
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time-dependent detection operators 

c s (t s ) = U\T)c s U(T). (44) 

Their action on the density matrix follows from eq. (21) when one uses that c a (t) corresponds 
to the direction ~it a (t) = —ysinSt + 2; cos 5t, and c&(£) to the opposite direction libit) = 
y sin 5t — 'z cos 5t. This gives 

c s (t s )p(R, 0, M(ts) = R 2 g s (6, 0, t s )p(R, 0, 0), (45) 

with g s (9, (j),t) = (1 + ~Tt ■ ~Tt s (t))/2. The general expression (17) for the contribution to the 
density matrix from a detection history {t iy s,} with the initial state (18), is found as 

L 

p L ({U, s t } , T) = exphiPa - e~ rT )] JJ(ri? e - rt ') 

x f d^f{9A)\[[gsXd^MU{T)me- TT/ \e^)U\T). (46) 
J 1=1 

Each detection s leads to a multiplication of the distribution function over the Bloch sphere 
by a factor g s {6, <fi, t) that now depends on the detection time. This time dependence corre- 
sponds to a rotation of the direction ~Tt s in the yz-plane. 

For the initial state of two decoupled modes, with a uniform distribution of the phase, 
the function / is uniform over the equator of the Bloch sphere. A detection at time t of a 
particle emitted by mode A or B then multiplies the distribution over the relative phase 
by the factor g a (4>) — (1 — sin5tsin0)/2 , or g b ((p) = (1 + sin<5tsin0)/2. These functions 
have their maximum value for = 3n/2 or = n/2. Strictly speaking, this distribution 
describes the state of the system in the Heisenberg picture, where it is not affected by 
continuous evolution, but only by the quantum jumps that describe the effect of detections. 
The evolution of the phase distribution during a typical detection history is conceptually 
simple. The total decay rate, summed over both detectors, is autonomous, and has the time 
dependent rate TR 2 exp(— Ft). The branching over the two detectors a and b is determined 
by the expectation value of g a (<f>) and gb(<f>), which has a contrast that oscillates in time at 
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the coupling frequency 5, as a result of the mode coupling. The effect of a detection on the 
phase distribution is a multiplication with the same factor (1 =p sin<5tsin0)/2, for detector 
a and b. This will eventually lead to convergence to the phase distribution to a single peak 
at a value where either one of the factors g s is maximal, hence = n/2 or = 3n/2. The 
convergence to these peaked distributions is slower than in the case of a detections through 
a single beam splitter, as a result of the oscillations of the contrast in the functions g s (t). In 
Fig. 4 we plot a set of typical phase distributions after L = 10 detections. The instants of 
detection are randomly selected, and the most probable dteection channel at that instant is 
chosen. The different curves correspond to a different selection of the instants of detection. 
As seen in Fig. 4, after each such history, the distribution over is a peak centered either 
at 7r/2 or at 37r/2. 



C. Coupling and energy shift 

An energy difference he between the two modes in addition to the effect of tunneling is 
described by the Hamiltonian 

H = -h5J x + heJ z , (47) 

which replaces (35). The angular-momentum operators are defined in eq. (1). We consider 
the same detection scheme used in the preceding subsection. The energy difference modifies 
the detection statistics and the phase distribution following a representative detection his- 
tory. On the Bloch sphere, the modified evolution operator U (t) is represented by a rotation 
in the positive direction around the axis e'z — Sx, over an angle Qt, with Q = y/e 2 + 5 2 . Equa- 
tions (43) for the density matrix after a detection history and (44) for the detection operators 
in the Heisenberg representation c s (t) remain valid. The detection operators are represented 
by points Tt s on the sphere that are reached from the poles when the opposite rotation is 
applied. Since the rotation axis does not lie in the equator plane, the azimuthal angle varies 
continuously with time, and the relative phase is no longer projected preferentially onto the 
same value. These unit vectors are found in the form 
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... , . , e5 . _ . ^ 5 . ^ / <5 2 _ \ ^ 
tf a (i) = - it 6 (i) = — (cosiM - l)x - —smilty + I — cosi« + — I z. 

They determine the factors g s (6,<j),t) = (1 + Tt ■ Tt s (t))/2 that multiply the distribution 
over the sphere when a particle emitted by mode A or B is detected. 

As above, we consider the case of an initially factorized state, which is represented by a 
uniform distribution over the equator of the Bloch sphere. When a particle from mode A or 
B is detected, the distribution over is multiplied by 

9a{4>) = ^ ^1 + ^ cos 0(cos fit - 1) - ^ sin sin VLt^j , 

1 ( sS $ 
9b{4>) = — I 1 — tz^ cos 0(cos Vtt — 1) + — sin sin Vtt 

The maximum of these functions no longer coincide with the maximum of ± sin 0, as is the 
case when e — 0. 

In Fig. 5 the resulting phase distributions are shown after a number of typical detection 
histories, each consisting of 10 detections, for e/S — 1/4. The prescription of the calculation 
is the same as used in Fig. 4. Now not only the width of the peak, but also their position 
varies for different selections of the detection times. This can be explained from the variation 
in the position where the maximum of g s (<j),t) occurs. 



V. LINEAR AND CIRCULAR CHAINS OF MODES 

The dynamics of a coupled chain of condensates in an optical lattice has been explored, 
with emphasis on the difference between a linear and a circular chain [16]. The coupling 
was due to tunneling between neighboring modes. One expect anaologous differences in 
the situation considered in this paper, where the phase relation between neighboring modes 
arises by spontaneous symmetry breaking from the observation of emitted bosons interfering 
through a beam splitters. This raises the question of the transitivity of the relative phase. 
When the relative phase between two modes A and B is well-determined, and the same holds 
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for the relative phase between two modes B and C, then one expects the phase between C 
and A should also be fixed. On the other hand, when this latter phase is also selected by 
direct interaction, one may expect different dynamics depending on whether the two paths of 
phase determination converge to the same result or not. In the present section we compare 
the phase dynamics on a linear and a circular chain of modes. 



A. Linear chain of modes 

We consider a linear chain of modes, as sketched in Fig. 6. As initial state we take the 
uncorrelated state given by the factorized density matrix 

p(0) = Yl Ps = ps-i ® Ps ® ps+i • • • , (48) 



s 



where the density matrix p s of each mode s has the form (34) with a uniform phase <j> 8 . 
Beam splitters are mixing the bosons emitted from neighboring modes s and s + 1, with 
orthogonal detection operators in the output channels 

d s ± = -j=(a a ± e-*'a a+1 ). (49) 

with Oj the annihilation operator of mode i. The evolution is described by the master 
equation (12), with 

r 

£>P= - S ^^&sp + pa\a s ),C\ = ^2(£is+ + Cis-), (50) 

s s 

where the contribution to C\ corresponding to the detection channels s± are specified by 

£is± = -d a ±pdl±. (51) 

Physically it is obvious that the detection statistics over the output channels of each beam 
splitter is identical to the statistics for each of the two beam splitters in Sec. Ill, since each 
mode emits into two input channels with equal rate. The density matrix corresponding to a 
given detection history with n s detections in channel s+, and m s detections in the channel s_ 
is easily written down by using that a detection in channel s + gives a factor cos 2 ((<3> s — £ s )/2), 

22 



and a detection in channel s_ a factor sin 2 (($ s — £ s )/2). After each detection history, the 
distribution over the phases <p s of all modes factorizes into a product of distributions for 
each relative phase $ s = (f) s — <p s+ i between neighbors. After n s detections in channel s+, 
and m s detections in the channel s_, the distribution over the relative phase <p s — s+ i is 
proportional to cos 2ns (($ s — £ a )/2) sin 2m-s ((<I> s — £ s )/2), and the distribution over the phases 
is proportional to the product 

n[cos 2 -(^)sin 2 -(^)]. (52) 

s 

Because of this factorization, the detection statistics for the pair of output channels of each 
beam splitter is uncorrelated to the other detections. The total number M s of detections in 
the time interval [0, T] on the two output channels of a single beam splitter is Poissonian 
with average value r 2 [l — exp(— TT)], and the probability distribution of the M s detections 
over the two detectors is identical to the distribution (33) [10]. Therefore, the most probable 
histories with M s detections on this sth beam splitter are given as {n si m s ) = (M s ,0) and 
(0, M s ). The relative phase <3> s between modes s and s + 1 converges to a single peak located 
at £ s or £ s + n, for each value of s. This also determines in a unique and unambiguous way 
the relative phase between any pair of modes. Hence for a linear chain of modes, the relative 
phase between two neighbors converges to one out of two possible values, in precisely the 
same way as it occurs for two modes and a single beam splitter. Spontaneous symmetry 
breaking occurs independently for each neighboring pair. 



B. Circular chain of modes 

Now we consider a series of K modes, coupled by beam splitters, and arranged into a 
circular chain. For K — 3, the scheme is presented in Fig. 7. Equations (48)-(50) still 
hold, with the index s running from 1 to K. The relative phases and the detection 
operators d s ± are defined as above for s — 1, 2,. . ., K — 1, while we denote $x = <pK — </>i, 
dx± = (flK ± e~^ K ai)/y/2. The number of beam splitters is now equal to the number of 
modes. On the other hand, since 
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K 

$>s = 0, (53) 

s=l 

the K modes have only K — 1 independent relative phases $ s , which makes the detection 
system overdetermined. This is the main difference with the case of the linear chain. Detec- 
tions on the sth beam splitter tend to drive the relative phase $ s to the value £ s or £ s + ir. 
However, these values are consistent only when the values of all £ s add up to a multiple of ir. 
The probability p({n s ,m s },T) of a specified number of detections by each detector in the 
time interval [0, T] factorizes as in eq. (28) in a Poisson distribution for the total number L 
of detections, with the mean value Kr 2 {l — e~ rT ), and the probability pl{{ti s , m s }) that the 
L detections are distributed over the detectors according to the indicated partition. This 
latter distribution can be specified in analogy to (29) by 

D 

P L ({n s , m s }) = wF({n s , m s }) (54) 



with 

F({n s ,m s })= \^—J J d0id0 2 JJ 



COS ( ) sin 3 { 



(55) 



After a detection history with n s detections in channel s+, and m s detections in the channel 
s_, the distribution over the relative phase is still proportional to (52). However, because of 
the relation (53), the relative phases are no longer independent, and the detection statistics 
of the output channels of the different beam splitters become correlated. 

The most probable histories can now be found by similar considerations as we used above 
in Sec. HID. For a total number of L = K x M detections, one might expect that the same 
number (M) of particles reaches each beam splitter, with the partition (n s ,m s ) = (M, 0) or 
(0, M) for all of the K beam splitters. This would indicate that the corresponding relative 
phases probed by these beam splitters will have converged to the value £ s or £ s = ir. However, 
in general this can only be true for all relative phases except one, because of the phase relation 
(53). Assume that this excepted relative phase has the index sq- As a result of this relation, 
the value of the last relative phase $ So is thereby also fixed. The distribution over the two 
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output channels sq+ and sq- will then be binomial, and the most probable partition is given 
by (n S0 ,m S0 ) = (M cos 2 (($ so - f so )/2), Msin 2 (($ so - £ so )/2). For symmetry reasons, each 
beam splitter has the same probability to end up in such a binomial distribution rather than 
a bunching one. The situation can be summarized by stating that in addition to the local 
spontaneous symmetry breaking for each beam splitter, also a global symmetry breaking 
occurs, by which the relative phase between two neighbors is not determined by the setting 
of their own shared beam splitter, but by the settings of all the other ones. 

We backed up this conclusion by a numerical calculation in the case of a three-mode 
ring, as sketched in Fig. 7. The settings of the beam splitters are given by £1 = £2 = 0, 
£3 = 7r/2. After 30 detections, one of the partitions with the highest probability was found 
to be (rii,mi) = (5,5), (n 2 ,m 2 ) = (10,0), (n 3 ,m 3 ) = (10,0). As one would expect from 
symmetry considerations, other partitions with the same maximal probability are found by 
swapping n s and m s for each beam splitter s, and also by a permutation of the three indices 
1, 2 and 3. 

VI. DISCUSSION AND CONCLUSIONS 

The absolute phase of a single-mode or multimode bosonic system is fully undetermined 
when the state of the system is diagonal in the total particle number. For bosonic atoms, 
this must be the case, since states with different particle numbers do not superpose. For 
a two-mode system we use the Schwinger representation with fictitious angular momentum 
operators to take advantage of the underlying SU(2) symmetry of the state space. This 
allows us to represent the density matrix of the two-mode system with an undetermined 
absolute phase and a Poissonian distribution of the total number of particles as an integral 
over the Bloch sphere of the fictitious angular momentum. The representation is given in 
eq. (10), where f(6,<f>) is the distribution function over the sphere. It may be viewed 
as the Glauber-Sudarshan P function restricted to the sphere. The azimuthal angle is 
the relative phase, whereas the polar angle 6 measures the ratio of the average number of 
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particles in A and B, with equal populations represented by points on the equator, and 
the poles representing states with all particles in one mode. The merit of these states with 
Poissonian distribution of the total particle number is that the overall decay of the modes 
factors out, and the detection statistics is the product of time-dependent probabilities for 
the total number of detections, and time-independent distributions for the partitions over 
the various detection channels. The effect of a detection is described by the action of an 
annihilation operator, which also corresponds to a point on the sphere. This is equivalent 
to the multiplication of the distribution function / (9, 0) by a factor that depends only on 
the distance over the sphere between the points (9, 4>) and the detection point. This allows 
exact expressions, both for the detection statistics, and for the conditional density matrix 
of the system for a given detection history. It also implies that identical detection statistics 
arises for different choices of the distribution / and the detection points on the Bloch sphere, 
provided that the setup has the same relative geometry on the sphere. This can correspond 
to quite different experimental setups, since the effect of detection through a beam splitter 
can be produced by a pulsed tunneling coupling between the modes. 

In the case that the modes are constantly coupled by tunneling, and in the presence of 
an energy difference between the modes, the phase distribution still becomes non-uniform 
by the detecting particles emitted by the two modes. However, since the preferred phase 
imposed by the detections is not the same for all detections in this case, the maximum in 
the phase distribution will continue to vary in position even after many detections. The 
convergence of the phase will be perturbed more strongly when interparticle interactions are 
important during a detection history. 

We treat explicitly the case of two modes which both emit particles in an input channel 
of two different beam splitters. When the settings of the beam splitters are different, they 
can drive the relative phase of the modes to values which are conflicting. Such a situation of 
conflicting phase values occurs for any number of modes which are coupled by beam splitters, 
and arranged in a circular chain. Our model shows that in these cases the most probable 
detection histories lead for each pair of neighboring modes to a relative phase converging 
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with equal probability to one of the conflicting values. The partition of the detection over 
the channels is a signature of the location of the peak in the phase distribution. Such a 
conflict does not arise for a linear chain of modes coupled by a beam splitter. A common 
feature of these various cases is that an initially factorized state of several modes builds up a 
specific value of all relative phases by only detecting their decay products in interference. In 
principle, this means that the modes become entangled, even though they have never been 
in direct contact. 
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FIGURES 




FIG. 1. Sketch of setup with two decaying modes, where each mode emits particles into the 
input port of two different beam splitters. 
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FIG. 2. The probability distribution pL({n s }) versus m and for equal number of detections 
through both beam splitters. Here L = 40 and £ = ir/2. The most probable histories are marked. 
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FIG. 3. (a) emitted particles are detected directly, without the use of beam splitters; (b) 
emitted particles re detected through a beam splitter. For each case, the position of the detectors 
on the Bloch sphere, and the distribution of the state before detection are also shown. 
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FIG. 4. Relative phase distributions for two coupled modes after L = 10 detections. Each curve 
of the ten curves corresponds to a different realization of the randomly selected detection times. 
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FIG. 5. Same as fig. 4, now for two coupled modes at different energy. Coupling strength and 
energy splitting specified by e/6 = 1/4. 
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FIG. 7. Sketch of setup with 3 modes arranged on circular chain. 
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